#### Set up ####
remove(list = ls())
library(tidyverse)
library(haven)
library(y2clerk)
library(extrafont)
#### Read in data ####
step <- read_dta("data/gss7221_r1a.dta")
responses <- step
responses[] <- lapply(step, unclass)
responses <- as.data.frame(responses)
rm(step)
# responses <- gss

#### Clean data ####
responses <- responses %>% 
  mutate(cappun_code = case_when(cappun == 2 ~ 0,
                                 cappun == 1 ~ 1))

#### Code RELTRAD Variable FROM ONLINE ####

library(car)
library(descr) #Get the CrossTable Function! Weighted! crosstab

#Get rid of the black oversamples - these throw off the proportions
responses = responses[responses$sample < 4| responses$sample==6| responses$sample>7,]
responses = as.data.frame(responses)
#recode into 5 major categories of religious affiliation
# 1) Protestant [Ask DENOM]	1371	47.8
# 2) Catholic
# 3) Jewish	
# 4) None	
# 5) Other (specify)
# 6) Buddhism
# 7) Hinduism
# 8) Other Eastern religion
# 9) Muslim/Islam
# 10) Orthodox Christian
# 11) Christian
# 12) Native American
# 13) Inter-/non-denominational
# 98) Don't know
# 99) No answer

responses$test=NA
responses$reltrad=NA
responses$xaffil = car::recode(responses$relig, "1=1;2=4;3=5;4=9;5:10=6;11=1;12=6;13=1")

responses$xaffil = as.factor(responses$xaffil)
levels(responses$xaffil) = c("prot", "cath", "jew", "other", "nonaf")

# The following code breaks down religious groups by evangelicals, black
# Protestants, mainline, liberal and conservative nontraditional,
# and Protestant nondenomination/no denomination.

#####Black Protestants
#Create a racial indicator
responses$black = ifelse(responses$race == 2, 1, 0)
responses$white = ifelse(responses$race == 1|responses$race == 3, 1, 0)
#Take the "other" Protestant denominations and pull out the 
#historical Black denominations, e.g. COGIC
responses$xbp = responses$other

responses$xbp = ifelse(responses$xbp %in% c(7, 14, 15, 21, 37, 38, 56, 78, 79, 85, 86, 
                                            87, 88, 98, 103, 104, 128, 133), 1, 0)


#National baptists and AME, AMEZ
responses$xbp = ifelse(responses$denom %in% c(12, 13, 20, 21), 1, responses$xbp)

#Blacks in certain denoms get recoded as Black Protestant
#Other baptist, amer. baptist, south. bap, other Methodists
# Create indicator for black race in denom and other
responses$bldenom = 0
responses$bldenom = responses$denom * responses$black
responses$blother = responses$other==93 * responses$black
#Call these black denominations
table(responses$bldenom[responses$year==2018])

responses$xbp = ifelse(responses$bldenom %in%
                         c(23, 28, 18, 15, 10, 11, 14), 1, responses$xbp)
#Black missionary baptists
responses$xbp[responses$blother == 1] = 1

#Evangelical Protestants#
#Recode the evangelicals in the other variable
responses$xev=responses$other
evother=c(2, 3, 5, 6, 9, 10, 12, 13, 16, 18, 20, 22, 24, 26, 27, 28, 31, 32,
          34, 35, 36, 39, 41, 42, 43, 45, 47,51, 52, 53, 55, 57, 63, 65, 66, 
          67, 68, 69, 76, 77, 83, 84, 90, 91, 92, 94, 97, 100, 101, 102, 106,
          107, 108, 109, 110, 111, 112, 115, 116, 117, 118, 120, 121, 122, 124, 
          125, 127, 129, 131, 132, 134, 135, 138, 139, 140, 146)
responses$xev=ifelse(responses$xev %in% evother,1,0)
#Cons Lutherans, cons presbyterians
responses$xev=ifelse(responses$denom %in% c(32,33,34,42), 1, responses$xev)
#White baptists, white other methodists
responses$whdenom = responses$denom * responses$white

responses$xev = ifelse(responses$whdenom %in%  c(10,18,15,23,14),1,responses$xev)
#Missionary baptist
responses$wother = responses$other * responses$white
responses$xev[responses$wother==93] = 1
responses$xev[responses] = 0
responses=as.data.frame(responses)
#Lifeway correction to reltrad
responses$xtn = responses$relig
responses$denom2 = responses$denom
#70 = No denomination or non-denominations
responses$denom2 = car::recode(responses$denom2, "70=1; else=0")
responses$xtn = car::recode(responses$xtn, "11=1; else=0")
responses$xtn[responses$denom2 == 1] = 2
responses$xtn = car::recode(responses$xtn, "1=1; 2=0")
#Only weekly or +weekly attenders
responses$xtn[responses$attend < 4|responses$attend==3|responses$attend==0|is.na(responses$attend)] <- 0
responses$xev[responses$xtn ==1] <- 1

responses$inter <- responses$relig
#Interdenominationals
responses$inter <- car::recode(responses$inter, "13=1; else=0")
responses$inter[responses$attend < 4|responses$attend==3|responses$attend==0|is.na(responses$attend)] <- 0
responses$xev[responses] <- 1

# Mainline Protestants
responses$xml = NA
responses$xml = responses$other
mpother=c(1,8,19,23,25,40, 44, 46, 48, 49, 50, 54, 70, 71, 72, 73, 81, 
          89, 96, 99, 105, 119, 148)
responses$xml=ifelse(responses$xml %in% mpother,1,0) 
#The denom category
responses$xml = ifelse(responses$denom %in% 
                         c(30, 50, 35, 31, 38, 40, 48, 43, 22, 41),1,responses$xml)

#Mainline baptist denom and methodists - if the R is white, they get coded mainline
responses$xml = ifelse(responses$whdenom %in% c(11, 28),1,responses$xml)

#Catholics
responses$xcath = responses$other 
#Polish National Church and Catholic
responses$xcath = ifelse(responses$other %in% c(123, 28),1,0)
#People who say that they are other get coded zero
responses$xcath=ifelse(responses$xaffil=="cath", 1, responses$xcath) 

#Jews
responses$xjew=0 
responses$xjew=ifelse(responses$xaffil=="jew",1,0)

#Adherents of other religions.
responses$xother = responses$other
responses$xother = ifelse(responses$xother %in% 
                            c(11, 17, 29, 30, 33, 58, 59, 60, 61, 62, 64, 74, 75, 80, 
                              82, 95, 113, 114, 130, 136, 141, 145),1,0)
#Adds others from main religious recoding
responses$xother=ifelse(responses$xaffil=="other" & responses$xev==0,1,0)

#Unaffiliateds/Nonaffiliateds
responses$xnonaff=0
responses$xnonaff[responses$xaffil=="nonaf"]=1

#The recodes non-denoms based on their attendance 
#Non active Don't Know Protestants coded to nonaffil
responses$xprotdk = ifelse(responses$denom == 70,1,0)
responses$xprotdk[responses$xprotdk == 1 & responses$attend >= 4] = 0
#Active Don't Know Protestants coded to evangelicals
responses$xnonaff[responses$xprotdk]=1
responses$xev[responses$xprotdk == 1 & responses$attend >= 4] = 1

#All these folks get coded evangelical
responses$reltrad = factor(NA, levels=c("Conservative Protestant",
                                        "Mainline Protestant",
                                        "Black Protestant",
                                        "Roman Catholic",
                                        "Jewish",
                                        "Other",
                                        "None"))

responses$reltrad[responses$xev==1]="Conservative Protestant"
responses$reltrad[responses$xml==1]="Mainline Protestant"
responses$reltrad[responses$xbp==1]="Black Protestant"
responses$reltrad[responses$xcath==1]="Roman Catholic"
responses$reltrad[responses$xjew==1]="Jewish"
responses$reltrad[responses$xother==1]="Other"
responses$reltrad[responses$xnonaff==1]="None"
responses$year = as.factor(responses$year)
responses = as.data.frame(responses)
# save(gss,file="gss7216_reltrad.csv")
#End of my poorly written R code! Sorry - I'll clean it up some day!

#### Relevel RELTRAD order ####
responses <- responses %>% 
  mutate(reltrad = fct_rev(reltrad))
#### Make plot ####
responses <- responses %>% 
  mutate(reltrad = fct_rev(reltrad))

relig_sums <- responses %>% 
  group_by(year,reltrad) %>% 
  summarize(Mean = mean(cappun_code, na.rm = T)) %>% 
  filter(!is.na(reltrad)) %>% 
  filter(reltrad != "Jewish",
         reltrad != "Other") %>% 
  filter(year != 1972) %>% 
  mutate(year = as.numeric(as.character(year)))

relig_sums$reltrad <- ordered(relig_sums$reltrad, levels = c("Conservative Protestant",
                                                              "Mainline Protestant",
                                                              "Roman Catholic",
                                                              "None",
                                                              "Black Protestant"))

ggplot(relig_sums, 
       aes(x = year,
           y = Mean)) +
  geom_line(aes(group = reltrad,
                linetype = reltrad)) +
  scale_y_continuous(labels = scales::label_percent(accuracy = 1L)) +
  scale_linetype_manual(values=c("solid", "twodash", "longdash", "dotted", "dotdash")) +
  scale_x_continuous(breaks = seq(1974, 2018, by = 2))  +
  labs(caption = "Source: General Social Survey") +
  xlab("") +
  ylab("Percent supporting capital punishment") +
  theme_classic() +
  theme(text = element_text(family = "Times New Roman"),
        plot.caption = element_text(hjust = 0),
        legend.text = element_text(size = 7),
        legend.position = "bottom",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text.x = element_text(angle = 315),
        axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0)),
        legend.title= element_blank())

ggsave("figures/cappun-support-byreltrad.pdf", device = "pdf", width = 6, height = 4, units = "in", dpi = 600)

#### Racial Resentment over time ####
responses <- responses %>%
  mutate(rr2_lackwill = case_when(racdif4 == 1 ~ 1,
                                  racdif4 == 2 ~ 0),
         rr3_discrimination = case_when(racdif1 == 1 ~ 0,
                                        racdif1 == 2 ~ 1))


responses <- responses %>% 
  mutate(rr_single_binary = case_when(wrkwayup <= 2 ~ 1,
                                      wrkwayup > 2 ~ 0)) %>% 
  mutate(rr_index_new = rr2_lackwill + rr3_discrimination + rr_single_binary)



hist(responses$rr_index_new)
mean(responses$rr_index_new, na.rm = T)

rr_sums <- responses %>% 
  group_by(year, reltrad) %>% 
  summarize(Mean = mean(rr_index_new, na.rm = T)) %>% 
  filter(!is.na(reltrad)) %>% 
  filter(reltrad != "Jewish",
         reltrad != "Other") %>% 
  mutate(year = as.numeric(as.character(year))) %>% 
  filter(year >= 1994)

rr_sums$reltrad <- ordered(rr_sums$reltrad, levels = c("Conservative Protestant",
                                                             "Mainline Protestant",
                                                             "Roman Catholic",
                                                             "None",
                                                             "Black Protestant"))



ggplot(rr_sums, 
       aes(x = year,
           y = Mean)) +
  geom_line(aes(group = reltrad,
                linetype = reltrad)) +
  scale_linetype_manual(values=c("solid", "twodash", "longdash", "dotted", "dotdash")) +
  scale_x_continuous(breaks = seq(1994, 2018, by = 2)) +
  scale_y_continuous(breaks = seq(0, 3, by = .25)) +
  labs(caption = "Source: General Social Survey") +
  xlab("") +
  ylab("Mean racial resentment scores") +
  theme_classic() +
  theme(text = element_text(family = "Times New Roman"),
        plot.caption = element_text(hjust = 0),
        legend.text = element_text(size = 7),
        legend.position = "bottom",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text.x = element_text(angle = 315),
        axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0)),
        legend.title= element_blank())

ggsave("figures/rr-byreltrad.pdf", device = "pdf", width = 6, height = 4, units = "in", dpi = 600)

        